//K K Holsman Sept 2011
//KIR'S HALF BAKED ADMB CODE
TOP_OF_MAIN_SECTION
  arrmblsize = 10000000;
DATA_SECTION
// !! ad_comm::change_datafile_name("M2.dat");
  init_int npred;  // number of predator spp
  init_int nprey;  // number of prey spp
  init_int nspp;   // number osrv_age_obsf predator spp
  init_int nyrs                // number of years
  init_ivector nages(1,nspp)   // number of unique legnth classes for each spp
  init_int nlengths
  init_ivector lengths(1,nlengths)          // legnths
  init_3darray N(1,nspp,1,nyrs,1,nlengths)          // number of obs in each  length matrix for each spp length classes vector
 int maxL
 int predd
 int preyy
 int j
 int k
 int i
 int aa
  !! maxL=nlengths;

 vector switch1(1,maxL) //pre-allocate a vector of 1s 
 vector preyuse(1,maxL)

 init_matrix Kglm_a(1,nspp,1,nspp)
 init_matrix Kglm_b(1,nspp,1,nspp)
 init_int nTyrs
 init_vector Tyrs(1,nTyrs)
 init_vector TempC(1,nTyrs)
 init_vector S_a(1,nspp)
 init_vector S_b(1,nspp)
 init_vector S_b2(1,nspp)
 init_vector aLim(1,nspp)
 init_vector bLim(1,nspp)
 init_vector aLW(1,nspp) 
 init_vector bLW(1,nspp)
 init_matrix maxK(1,nspp,1,nspp)
// init_3darray N(1,nspp,1,nyrs,1,maxL)
// init_imatrix AvgN(1,nspp,1,maxL)
 init_matrix LA_age(1,nspp,1,nages)// length to age coversion matrix
 init_matrix LA_Lengths(1,nspp,1,nages)// length to age coversion matrix

  
  matrix AvgN(1,nspp,1,maxL)  //pre-allocate mean Number per length
  matrix S(1,npred,1,maxL)  //pre-allocate mean stomach weight as a function of age
  matrix Limit(1,npred,1,maxL)  //pre-allocate
  3darray K(1,npred,1,nprey,1,maxL) //K(1) is npreyXmaxL matrix of stomach proportions of predator 1 
  //Limit<-matrix(0,npred,nlengths[1])//pre-allocate
  4darray U(1,npred,1,nprey,1,maxL,1,maxL) //pre-allocate U(pred#,Prey#,predsize,preysize), so that U(2,1)= nlengths(2)x nlengths(1) matrix 
  //U(1,pred,1,nprey,1,predAge,1,preyAge)
  //representing predator 2, and prey 1. The proportion of each prey size in the diet of each predator size.
  3darray ration(1,npred,1,nyrs,1,maxL)
  3darray Consum(1,nyrs,1,npred,1,maxL) //pre-allocate, indiviudal consumption in grams per predator
  3darray Esum(1,nprey,1,nyrs,1,maxL) //pre-allocate total grams of that prey eaten in a given year
  3darray    E(1,nprey,1,nyrs,1,maxL) //pre-allocate total grams of a prey eaten by each predator in a given year y, E(1) is a nrys x maxL matrix of grams eaten of prey 1
  3darray Eage(1,nspp,1,nyrs,1,nages)
  3darray Bage(1,nspp,1,nyrs,1,nages)
  3darray M2(1,nspp,1,nyrs,1,nages) //pre-allocate      total predation mortality in a given year for each prey, M2(1) is a nrys x maxL matrix of mortality for prey 1
  //matrix Nuse(1,npred,1,maxL) // pre-allocate the N of preds used in consumption
  3darray Nuse(1,nyrs,1,nspp,1,maxL)// pre-allocate the N of preds used in consumption, where Nuse varies by year
  matrix LA_Lengthshat(1,nspp,1,nages);// pre-allocate the N of preds used in consumption, where Nuse varies by year

PARAMETER_SECTION
  // need to convert maxL to nlengths to make into a ragged array - ask jim
//  init_bounded_matrix LA_Lengthshat(1,nspp,1,nages,0,200,2);
  init_vector a(1,nspp,1);
  init_vector b(1,nspp,1);
  objective_function_value f;

PRELIMINARY_CALCS_SECTION
////////////////////////////////////////////////////////////////////////////////////////////////////////
//// Calculate S: mean stomach weight (as prop of body weight) as a function of length 
////////////////////////////////////////////////////////////////////////////////////////////////////////

  for (predd=1;predd<=npred;predd++)
  {
    AvgN(predd)=N(predd,1);
    //  S[predd,]<-S_a[predd]+S_b[predd]*age_size[predd,]+S_b2[predd]*(age_size[predd,]2)
    S(predd)=S_a(predd)+S_b(predd)*LA_Lengths(predd)+S_b2(predd)*pow(LA_Lengths(predd),2);
    for (j=1;j<=nlengths;j++)
    {    
      if(LA_Lengths(predd,j)>80)
      {
        //S[predd,age_size[predd,]>80]<-S_a[predd]+S_b[predd]*80+S_b2[predd]*(80^2)
        S(predd,j)=S_a(predd)+S_b(predd)*80+S_b2(predd)*pow(80,2); // set everything above 80 to 80
      }
    }
  } 
INITIALIZATION_SECTION

PROCEDURE_SECTION
  cal_Nuse();
  calc_K();
  calc_lim();
  calc_U();
  calc_Consum();
  calc_E();
  calc_M2();
  dummy_fun();
  for (predd=1;predd<=npred;predd++)
  {
      f+= norm2(LA_Lengthshat(predd)-LA_Lengths(predd));
  }
  if (mceval_phase()) mcOutput();

FUNCTION dummy_fun
  for (predd=1;predd<=npred;predd++)
  {
    for (j=1;j<=nages(predd);j++)
    {
      LA_Lengthshat(predd,j)=value(exp(a(predd)+b(predd)*log(LA_age(predd,j))));
    }
  }
FUNCTION cal_Nuse
// recalculate Nuse for consumption with new estimates of N of each pred:
  for (i=1;i<=nyrs;i++)
  {
    for (predd=1;predd<=npred;predd++)
    {
      for (aa=1;aa<=maxL;aa++)
      {
        Nuse(i,predd,aa)=0;
        if(aa<=nlengths)
        {
          Nuse(i,predd,aa)=N(i,predd,aa);
        }
      }
    }
  }
// Nuse where a > nlengths(pred) are just set to 0
////////////////////////////////////////////////////////////////////////////////////////////////////////
// Calculate K (diet pref) for each size of predator - can be changed from sizes to ages using median or mean size per age
////////////////////////////////////////////////////////////////////////////////////////////////////////
FUNCTION calc_K
  for(predd=1;predd<=npred;predd++)
  {
    dmatrix Ktmp(1,nprey,1,nlengths);////pre-allocate matrix of 0 - need to fix
    Ktmp=1.0; // assigns value 0 to Ktmp
    for(preyy=1;preyy<=nprey;preyy++)
    {
      for(j=1;j<=nlengths;j++)
      {
        //Ktmp[preyy,j]<-min(1,exp(Kglm_a[predd,preyy]+(age_size[predd,j]*Kglm_b[predd,preyy])))
        Ktmp(preyy,j)=min(1.,exp(Kglm_a(predd,preyy)+(LA_Lengths(predd,j)*Kglm_b(predd,preyy))));//Ktmp is a matrix
        if(Ktmp(preyy,j)>maxK(predd,preyy))
        {
          Ktmp(preyy,j)=maxK(predd,preyy); // broken stick upper limit (see graph)
        }
        if(predd==1)
        {
          //if(predd==1){Ktmp[preyy,age_size[predd,]>90]<-0}
          if(LA_Lengths(preyy,j)>90)
          {
            Ktmp(preyy,j)=0;  // force the preference to zero (aka other prey) for pollock greater than 90 (see graph)
          }
        }
      } // end for each size/age                 
    } // end for each prey  
    K(predd)=Ktmp; // assign K of the array pred to the Ktmp matrix, repeat
  }//end pred
//        for(predd=1;predd<=npred;predd++)
//        {
//
//            for(j=1;j<=nlengths;j++)
//            {
//                dvariable Ktmp;
//                Ktmp=0;
//                for (k=1;k<=nprey;k++)
//                {
//                     Ktmp+=K(predd,k,j); // column sum for K
//                }// end for k nprey
//    if(Ktmp>1)
//                {
//                       //Ktmp[,colSums(Ktmp)>1]<-Ktmp[,colSums(Ktmp)>1]/colSums(Ktmp[,colSums(Ktmp)>1])  //re-standardize for columns where the sum is greater than 1
//          K(predd,k,j)= K(predd,k,j)/Ktmp;  //re-standardize for those  //re-standardize for those values that exceed 1 the vector Ktmp(j) is standardized by the rowsums
//    }// end if sum  
//            }// end for j ages
//        } // end for predd
                //K[,,predd]<-Ktmp

////////////////////////////////////////////////////////////////////////////////////////////////////////
// Calculate lim upper size limit that each predator can eat (not prey dependent) - can be changed from sizes to ages using median or mean size per age
////////////////////////////////////////////////////////////////////////////////////////////////////////

FUNCTION calc_lim
  for(predd=1;predd<=npred;predd++)
  {
   // upper limit of prey size that each pred of each size can consume
    Limit(predd) = aLim(predd)+bLim(predd)*lengths(predd);
  }//end pred

////////////////////////////////////////////////////////////////////////////////////////////////////////
////based on AGE Calculate (U) effective prey size freq (Npa/sum(N))*switch -- rows=prey, cols=preds, can be in calcs section if Npa/Np doesn't change each year
////////////////////////////////////////////////////////////////////////////////////////////////////////
FUNCTION calc_U
//U(1,pred,1,nprey,1,predAge,1,preyAge)
  for(predd=1;predd<=npred;predd++)
  {
    for(preyy=1;preyy<=nprey;preyy++)
    {
      dmatrix Utmp(1,nlengths,1,nlengths);
      Utmp=0; // pre-allocate matrix of 0, fix
      dvector preyprop(1,nlengths);
      preyprop=0;
                                               
      preyprop=AvgN(preyy)/sum(AvgN(preyy));////////// <- change this to numerical abundance by year !!!!  ? need to pre-allocate?
      for(j=1;j<=nlengths;j++)
      {// for each j size of preyy
        switch1=1; //pre-allocate a vector of 1s 
          preyuse=0;
        for(k=1;k<=nlengths;k++)
        {// for each k size of predd
          if(lengths(j)>=Limit(predd,k))
          {//if the prey size j is greater than or equal to the predator gape limit of pred size k....
            switch1(k)<-exp(-(lengths(k)-Limit(predd,k))/4);  // exp decline in preferance if above size limit
          }
        }//end k
                                
                                // ask jim: does the above equation need anything special to mulitiply two identical vectors?
        Utmp(j)=preyprop(j)*switch1*K(predd,preyy);
                                 //Utmp is a vector of 1xlengths (k) of pred lengths, switch is a 1xk vector of uses or not, and K is a 1xk vector of diet proportions
      }
      U(predd,preyy)=Utmp;// U tmp is a matrix of preysize X predsize U (npred,nprey,preysize,predsize)
    }//end prey
  }// end pred
  
////////////////////////////////////////////////////////////////////////////////////////////////////////
//// Calculate annual ration C=24*0.0134*exp(0.0115*TempC)*91.25*S  units are g of food/g of predator
////////////////////////////////////////////////////////////////////////////////////////////////////////
FUNCTION calc_Consum
  for(i=1;1<=nyrs;i++)
  {
    for(predd=1;predd<=npred;predd++)
    {
      for(j=1;j<=nlengths;j++)
                        {
           Consum(predd,i,j)=24*0.0134*mfexp(0.0115*TempC(i))*91.25*S(predd,j);
                             // vector of specific consumption rates for each size of predator in terms of g/g/yr of all prey consumed per predator per year
                 ration(predd,i,j)=Consum(predd,i,j)*Nuse(i,predd,j)*aLW(predd)*pow(j,bLW(predd));//*mnwt_bin(predd,j);
                              // need to change this to use the N from the previous year!
                             // annual ration for each size class of predator each year ration (pred,yrs,predages) ration(1) is nyrs x nlengths matrix for predator 1
                        }
                }
  }
////////////////////////////////////////////////////////////////////////////////////////////////////////
//// Calculate grams of prey eaten (E) units are g of prey E= U*C*N(pred,pred_size)*wt(pred,pred_size)
////////////////////////////////////////////////////////////////////////////////////////////////////////
FUNCTION calc_E
         //U(1,pred,1,nprey,1,predAge,1,preyAge)
         //E(1,npreyy,1,nyrs,1,preyAges)  //gives the number of each prey size eaten by all predators in the model for each year
 for(i=1;i<=nyrs;i++)
 {
   for(preyy=1;preyy<=nprey;preyy++)
   {  
     for(j=1;j<=nlengths;j++)
     {//for each prey length j
               dvector Etmp(1,nlengths);
               Etmp=0;
               for (predd=1;predd<=npred;predd++)
               {
                   Etmp(j)+=U(predd,preyy,j)*ration(predd,i);// Etmp(j) is a 1xk vector of biomass consumed for each pred size k
                                                                //U(predd,preyy,j) is a 1xk vector of suit. for each pred size k, 
                                                                //ration is a 1xk vector of pred rations
               }//end for predd
        E(preyy,i,j)=sum(Etmp);// Etmp is a vector of biomass consumed for all preds sizes 1-k, sum is the total biomass consumed for prey size j, in year i
        }// end for each prey size j
     }// end for prey preyy
  }//end for year i


////////////////////////////////////////////////////////////////////////////////////////////////////////
//// LASTLY!! Calculate M2 as E/N for each AGE
////////////////////////////////////////////////////////////////////////////////////////////////////////
FUNCTION calc_M2
  for(i=1;i<=nyrs;i++)
  {
    for(preyy=1;preyy<=nprey;preyy++)
    {
      Eage(preyy,i,1) = E(preyy,i,1);
      for(j=2;j<=nages(preyy);j++)
      {// for each age class of prey
        dvariable Eagetmp;
        Eagetmp=0;
        dvariable Bagetmp;
        Bagetmp=0;
        for(k=LA_Lengths(preyy,j-1)+1;k<=(LA_Lengths(preyy,j));k++)
        {// for each size k in the age class j of the prey 
          Eagetmp+=E(preyy,i,k);
          Bagetmp+=Nuse(i,preyy,k)*aLW(preyy)*pow(k,bLW(preyy));
        }
        Eage(preyy,i,j)=value(Eagetmp);  // Eage is the sum of the grams eaten for each of the sizes in the age group j
        Bage(preyy,i,j)=value(Bagetmp);  // Bage is the sum of the biomass for each of the sizes in the age group j
        M2(preyy,i,j)=Eage(preyy,i,j)/(Bage(preyy,i,j));  // need to make biomass updated in each year....
      }
    }
  }

////////////////////////////////////////////////////////////////////////////////////////////////////////
//// LASTLY!! Calculate M2 as E/N for each AGE
////////////////////////////////////////////////////////////////////////////////////////////////////////
//FUNCTION calc_M2
// for(i=1;i<=nyrs;i++)
// {
//     for(preyy=1;preyy<=nprey;preyy++)
//     {
//        //for(j=1;j<=nlengths(preyy);j++)
//        //{
//          //    biom(preyy,i,j)=Nuse(i,preyy,j)*mnwt_bin(preyy,j);
//        //}
//        
//  M2(preyy,i)=E(preyy,i)/(Nuse(i,preyy)*mnwt_bin(preyy));  // need to make biomass updated in each year....
//      }
// }
FUNCTION mcOutput
// Output parameters sampled; printed to 'cout' during the '-eval' phase.
// Pipe to a file, such as vbmc.dat. Each line has 6 values.
  cout << E << " " << ration << " " << S << " " << K << " " << U << " " << M2 << " " << Consum << endl;
REPORT_SECTION
   report << "E" << endl;
   report << E << endl;
   report << "Ration" << endl;
   report << ration << endl;
   report << "S" << endl;
   report << S << endl;
   report << "K" << endl;
   report << K << endl;
   report << "U" << endl;
   report << U << endl;
   report << "M2" << endl;
   report << M2 << endl;
   report << "Consum" << endl;
   report << Consum << endl;
   report << "$mcnames" << endl;
   report << "E ration S K U M2 Consum"   << endl;
   report << "$mcest" << endl;
   report << E << " " << ration << " " << S << " " << K << " " << U << " " << M2 << " " << Consum << endl;
